#include "ParticleSource1D_Emit.h"
#include "../Parallel/ParallelMPI1D.h"
#include "../Field/Field1D.h"

#include <cmath>
#include <iomanip>
#include <algorithm>
#include <ostream>
#include <sstream>

using namespace std;


ParticleSource1D_Emit::ParticleSource1D_Emit(PicParams* params_, ParallelMPI* pmpi_, PhysicalField* fields_, vector<Species*>* vecSpecies_, ParticleSourceParameter ps_param_):
    ParticleSource1D(params_, pmpi_, fields_, vecSpecies_, ps_param_)
{
    dt_ov_dx    = params->dt / params->cell_length[0];
    dt          = params->dt;
    YZArea      = 1.0;

    ySqrt_factor    = pow(params->const_e, 3.0) / (4.0 * params->const_pi * params->const_ephi0 * work_function *work_function);
    a_factor        = a_FN * params->const_e * params->const_e / work_function;
    b_factor        = -b_FN * pow(work_function, 1.5) / params->const_e;

    emit_step    = 1.0 + ps_param_.emit_particle_number * params->species_param[ps_param_.species_number].weight * params->cell_length[0] / ( ps_param_.emit_flux * params->dt );
    emit_rem     = ps_param_.emit_flux * emit_step * params->dt / ( params->species_param[ps_param_.species_number].weight * params->cell_length[0] ) - ps_param_.emit_particle_number;
    emit_rem_tot  = 0.0;
    log_info<<"emit_step:    "<<emit_step;
    log_info<<"emit_species: "<<ps_param.species_name<<" "<<ps_param.species_number;
}

ParticleSource1D_Emit::~ParticleSource1D_Emit()
{

}


void ParticleSource1D_Emit::emitLoad(int itime)
{
    Species1D   *species1D;

    if(itime%emit_step == 0)
    {
        double emit_number_temp;
        emit_rem_tot += emit_rem;
        emit_number_temp = ps_param.emit_particle_number;
        if(emit_rem_tot > 1.0)
        {
            emit_number_temp += 1;
            emit_rem_tot -= 1.0;
        }


        species1D = static_cast<Species1D*>((*vecSpecies)[ps_param.species_number]);

        ParallelMPI1D* pmpi1D = static_cast<ParallelMPI1D*>(pmpi);

        if(ps_param.emit_type == "fieldEmit")
        {
            //field emission is calculated using Fowler-Nordheim formulae
            //from "modelling vacuum arcs: from plasma initiation to surface interactions"
            Field1D<double>* Ex1D = static_cast<Field1D<double>*>(fields->Ex_);
            emitField = (*Ex1D)(0);
            //emitField *= params->norm_efield;
            emitJ = (a_factor * emitField * emitField / t_y2(emitField)) / exp(b_factor * v_y(emitField) / emitField);
            //emitJ /= params->norm_j;
            number_particle_emit = emitJ * dt_ov_dx / (weight_const * species1D->species_param.charge);
        }
        else if(ps_param.emit_type == "relEmit")
        {
            ParticleSource1D_Emit *relParticleSource_Emit = static_cast<ParticleSource1D_Emit*>(relParticleSource);
            number_particle_emit = relParticleSource_Emit->number_particle_emit * relEmit_factor;
        }
        else if(ps_param.emit_type == "regular")
        {
            number_particle_emit = emit_number_temp;
            log_debug<<"Injected number: "<<number_particle_emit;
        }

        if(ps_param.emit_position =="left" && pmpi1D->isWestern())
        {
            auto cell_particles = species1D->particles_[0];
            int number_particles_pre = cell_particles->size();
            cell_particles->add_particles(number_particle_emit);
            for(int i_particle = number_particles_pre; i_particle < cell_particles->size(); i_particle++)
            {
                auto& particle = cell_particles->data_[i_particle];
                particle.position[0] = (((double)rand() / RAND_MAX))*params->cell_length[0] * ps_param.emit_offset;

                double ran;
                do {
                    ran = (double)rand() / RAND_MAX;
                }
                while (ran == 0.0);
                //Velocity magnitude: from Maxwell velocity distribution
                //The angle between velocity of emitted particle and the surface normal: cosine
                //     cos(alpha) = sqrt(random number 0-1)
                //The azimuthal angle is uniformly distributed on the interval [0 2pi]
                double psm = sqrt(2.0 * params->const_e * ps_param.temperature / species1D->species_param.mass) * sqrt(-log(ran));
                double cosAlpha = sqrt((double)rand() / RAND_MAX);
                double sinAlpha = sqrt(1.0 - cosAlpha * cosAlpha);
                double phi   = 2.0 * M_PI*(double)rand() / RAND_MAX;

                particle.velocity[0] = abs( psm * cosAlpha );
                particle.velocity[1] = psm * sinAlpha * cos(phi);
                particle.velocity[2] = psm * sinAlpha * sin(phi);
            }
        }
        else if(ps_param.emit_position=="right" && pmpi1D->isEastern())
        {
            auto cell_particles = species1D->particles_[0];
            int number_particles_pre = cell_particles->size();
            cell_particles->add_particles(number_particle_emit);
            for(int i_particle = number_particles_pre; i_particle < cell_particles->size(); i_particle++)
            {
                auto& particle = cell_particles->data_[i_particle];
                particle.position[0] = params->sim_length[0] - (((double)rand() / RAND_MAX))*params->cell_length[0]*emitOffset;

                double ran;
                do {
                    ran = (double)rand() / RAND_MAX;
                }
                while (ran == 0.0);
                //Velocity magnitude: from Maxwell velocity distribution
                //The angle between velocity of emitted particle and the surface normal: cosine
                //     cos(alpha) = sqrt(random number 0-1)
                //The azimuthal angle is uniformly distributed on the interval [0 2pi]
                double psm = sqrt(2.0 * params->const_e * ps_param.temperature / species1D->species_param.mass) * sqrt(-log(ran));
                double cosAlpha = sqrt((double)rand() / RAND_MAX);
                double sinAlpha = sqrt(1.0 - cosAlpha * cosAlpha);
                double phi = 2.0 * M_PI*(double)rand() / RAND_MAX;

                particle.velocity[0] = -abs( psm * cosAlpha );
                particle.velocity[1] = psm * sinAlpha * cos(phi);
                particle.velocity[2] = psm * sinAlpha * sin(phi);
            }
        }
    }
}
